Lattice Melting and Rotation in Perpetually Pulsating Equilibria 
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Systems whose potential energies consists of pieces that scale as r~ 2 together with pieces that scale 
as r 2 , show no violent relaxation to Virial equilibrium but may pulsate at considerable amplitude for 
ever. Despite this pulsation these systems form lattices when the non-pulsational "energy" is low, 
and these disintegrate as that energy is increased. The "specific heats" show the expected halving 
as the "solid" is gradually replaced by the "fluid" of independent particles. The forms of the lattices 
are described here for N < 20 and they become hexagonal close packed for large N. In the larger 
N limit, a shell structure is formed. Their large N behaviour is analogous to a 7 = 5/3 poly tropic 
fluid with a quasi-gravity such that every element of fluid attracts every other in proportion to 
their separation. For such a fluid, we study the "rotating pulsating equilibria" and their relaxation 
back to uniform but pulsating rotation. We also compare the rotating pulsating fluid to its discrete 
counter part, and study the rate at which the rotating crystal redistributes angular momentum and 
mixes as a function of extra heat content. 
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INTRODUCTION 



Astrophysics has provided several new insights into 
ways statistical mechanics may be extended to cover a 
wider range of phenomena. Negative heat capacity, bod- 
ies which get cooler when you heat them, were first en- 
countered by [6] and when such bodies were treated ther- 
modynamically by [1][13], what seemed natural to as- 
tronomers was seen as an apparent contradiction in basic 
physics to those from statistical mechanics background. 
Even after a physicist [16] first resolved this paradox and 
[8] gave an easily soluble example, there was consider- 
able reluctance to accept the idea that micro canonical 
ensembles could give such different results to canonical 
ones. There was still some reluctance even in 1999 [9], 
though to those working on simulations of small clusters 
of atoms or molecules, the distinction was well under- 
stood by the early 1990s, see e.g. [3]. However by 2005, 
the broader statistical mechanics community had em- 
braced these ideas, and emphasized this difference which 
had been with us for 30 years (see Pichon & Lynden-Bell 
2006, where an early account of this work is given). An- 
other area to which statistical mechanics might extend 
is that of collisionless systems which may be treated as 
"Vlasov" fluids in phase space. While early attempts 
at finding such equilibria gave interesting formulae rem- 
iniscent of a Fermi-Dirac distribution [7], there is ample 
evidence from astronomical simulations that these equi- 
libria are not reached in gravitationnal systems. There 
are also different ways of doing the counting that lead to 
different results and recently, [2] gave an example that 
demonstrates that the concept of a unique final state de- 
termined by a few constants of motion found from the 
initial state is not realized. Thus the statistical mechan- 



ics of fluids in phase space remains poorly understood 
with theory and at best losely correlated with simula- 
tions and experiments (but see [4]). 

This paper is concerned with a third area of non stan- 
dard statistical mechanics which is restricted to very spe- 
cial systems, those for which the oscillation of the scale of 
the system separates off dynamically from the behaviour 
of the rescaled variables that are now scale free. These 
systems were found as a bi product of a study which 
generalized Newton's soluble N-body problem [10] [11] 
(papers I and II). A one dimensional model with exact 
solutions was found by [5] and the model considered here 
is a three dimensional generalisation of his combined with 
Newton's. A first skirmish with the statistical dynamics 
of the scale free variables despite the continuing oscil- 
lation of the scale was given there. Since then, [12], 
hereafter paper III, have showed that the peculiar ve- 
locities of the particles do indeed relax, as predicted, 
to a Maxwellian distribution, whose temperature con- 
tinuously changes as (scale) -2 . This occurs, whatever 
the ratio of the relaxation time to the pulsation period. 
The peculiar velocities are larger whenever the system is 
smaller. When the system rotates at fixed angular mo- 
mentum, the resulting statistical mechanics leads again 
to Maxwell's distribution function relative to the rotat- 
ing and expanding frame (see paper III, equation 17). 
It is then the distribution of the peculiar velocities after 
the "Hubble expansion velocity" and the time dependent 
rotation are removed that are distributed Maxwellianly; 
when the system is fluid, the density distribution is a 
Gaussian flattened according to the rotation that results 
from the fixed angular momentum. 

The interest of this problem for those versed in statis- 
tical mechanics is that it is no longer an energy that is 
shared between the different components of the motion. 
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5 Two tetrahedra joined on a face 

6 Two pyramids joined at their bases 

7 Two pentagonal pyramids joined at their bases 

8 A twisted cube 

9 Three pyramids, each pair sharing the edge of a base 

10 Skewed pyramid twisted 

11 Two skewed pyramid twisted and one center 

12 Icosahedron i.e. two skewed pentagonal pyramids 

13 Icosahedron and one center 

14 Hexagonal and pentagonal pyramids and one center 

15 Two skewed hexagonal pyramids and one center 

16 Pentagonal pyramid and an hexagon and a triangle and one center 

17 Hexagonal pyramid over an hexagon over a triangle and one center 

18 Two pentagonal pyramids and one twisted pentagon and one center 

TABLE I: First 18 configurations of equilibria. A few of them are shown on Figure 1. Note that as the number of elements 
increases, the final configurations need not be unique. 



The interest for astronomers lies in part because these 
systems suffer no violent relaxation to a size that obeys 
the Virial theorem. Nevertheless, despite the continual 
pulsation of such systems the rescaled variables within 
them do relax to a definite equilibrium. 

The form of interaction potential ensures that neither 
the divergence of the potential energy at small separa- 
tions nor the divergence of accessible volume at large 
separations which so plague systems with normal grav- 
ity occur for this model. Thus its interesting different 
statistical mechanics is simpler and free of any contro- 
versial divergences thanks to the small range repulsion 
and the harmonic long range attraction which extends to 
arbitrarily large separations. 

Although the harmonic long range attraction is not re- 
alised in nature (with the possible exception of quarks), 
nevertheless, within homogeneous bodies of elipsoidal 
shape ordinary inverse square gravity does lead to har- 
monic forces analogous to those found here. Section 2.1.1 
shows that only force laws with our particular scalings 
give such exact results. 

In this paper, we demonstrate that such systems can 
form solids (albeit ones that pulsate in scale). We study 
in Section III their behaviour in the large N limits and 
show that they stratify into shells. 

We show the phase transition as these structures melt 
and the corresponding changes in "specific heat". We 
also discuss the analogous a 7 = 5/3 fluid system (cor- 
responding to a classical white dwarf with an odd grav- 
ity, see below) , we predict their equilibrium configuration 
and investigate their properties when given some angular 
momentum. Section II derives the basic formulae for the 
N-body system and its fluid analog. 



II. DERIVATION 

A. The discrete system 

We consider N particles of masses rrii at , i = N and 
set 

Y j m i = M, ^m iri = Mr, (1) 

and half the trace of the intertial tensor as 

I = ^2m i (r i -r) 2 = Ma 2 , (2) 

which defines the scale a. In earlier work we showed 
both classically (paper I) and quantum mechanically [11] 
(paper II) that if the potential energy of the whole system 
was of the form 

V = W{a)+a- 2 W 2 (k), (3) 

where a is the 3N dimensional unit vector 
TV-i/2 

a = (n - r, r 2 - r .... , r N -r), (4) 

a 

then the motion of the scaling variable a separates dy- 
namically from the motions of both r and the r so those 
motions decouple. If we ask that V be made up of a sum 
of pairwise interactions then the hyperspherical poten- 
tial, W, (independent of direction in the 3N — 3 dimen- 
sional space) has to be of the form: 

W(a) = \?Ma 2 = V-2 , (5) 
where uo 2 is constant but depends linearly on M, and 

W 2 (a) =a 2 V 2 =a 2 ^^K 2 |r,-r J |- 2 . (6) 

i<j 

Hence W2 is a function of a and is independent of the 
scale a. 



3 




FIG. 1: (Color online) A small set of remarkable equilibria for low N; from left to right and top to bottom: N = [5] two 
tetrahedra joined at their bases, [6] two pyramids joined at their bases, [8] a twisted cube (three pyramids each pair sharing 
a base edge), [12] Icosahedron i.e. two skewed pentagonal pyramid, [13] Icosahedron+center, [14] hexagonal+ pentagonal 
pyramids+ center, [17] hexagonal pyramid over hexagon over triangle+center, [18] two pentagonal pyramids+twisted pentagon 
center. Java animations describing theses oscillating crystals are found at http://www.iap.fr/users/pichon/nbody.html. 



1. Relationship to the Virial Theorem 

Let us now see why potentials of the for Vz + V-2 are 
so special by looking at the Virial Theorem and the con- 
dition of Energy conservation. Take the more general 
potential energy to be a sum of pieces V n where each V n 
scales as r~ n on a uniform expansion, n can be positive 
or negative. Then V = ^2V n . For such a system the 
Virial Theorem reads: 

Now we already saw that V-2 oc I and V2 clearly drops 
out of the final sum because n — 2 is zero. Thus for 



potentials of the form V = V-2 + V2 

X -i =2E- W-2 = 2E- Alu 2 Q/^ , (7) 

so / vibrates harmonically with angular frequency 2uo 
about a mean value E/uj 2 and shows no Violent Relax- 
ation (1). Multiplying by 4/ and integrating 

P = 8EI - 4co 2 I 2 - 4M 2 C 2 , 

where the last term is the constant of integration chosen 
in conformity with Eq. (8). Now recall from Eq. (2) that 
/ = Ma 2 so we find on division by SM 2 a 2 that 

cl 2 jC 2 1 9 9 E , . 

Y + 2^ + r a =M> (8) 




are 



FIG. 2: (Color online) An example of large N static spher- 
ical equilibrium obeying Eq. (9); here only a given shell is 
represented for clarity. 



which we recognise as the specific energy of a particle 
with specific "angular momentum" C moving in a simple 
harmonic spherical potential of 'frequency' uo. In such a 
potential a vibrates about its mean at 'frequency' 2u. 

Note here that the 7 = 5/3 fluid also has a term 3(7 — 
1)U = 2U in the virial theorem since its internal energy, 
U scales like r~ 2 so it can be likewise absorbed into 
the total energy term. Thus if every elementary mass of 
such a fluid attracted every other with a force linearly 
proportional to their separation, then that system too 
would pulsate eternally, as described above in Eq. (7) 
(see Section II B below). 

The aim of this paper is to demonstrate the existence 
of perpetually pulsating equilibrium lattices and to study 
the changes as the non-pulsational 'energy' MC 2 /2 is in- 
creased. We show that the part of the potential left in 
the equation of motion for the rescaled variables is the 
purely repulsive V 2 . It is then not surprising that the lat- 
tice disruption at higher non-pulsational energy occurs 
quite smoothly and the solid phase appears to give way 
to the "gaseous" phase of half the specific heat without 
the appearance of a liquid with another phase transition. 
The hard sphere solid has been well studied and behaves 
rather similarly. 



2. The Equations of Motion and their separation 

Although the work of this section can be carried out 
when the masses ra$ are different, (see paper I), we here 
save writing by taking m % = m and so M = Nm. We 
start with V of the form (3). The equations of motion 



mvi — —dVjdvi 



(9) 



Now V is a mutual potential energy involving only — r j 
so ^dV/dri = and summing the above equation for 



all i we have 



d 2 r /dt 2 = 



r = r + ut 



Henceforth we shall remove the centre of mass motion 
and fix the centre of mass at the origin so r = 0. Now 
the 3N- vector a can be rewritten in terms of its length a 
and its direction a = a/a (a unit vector). Equation (9) 
takes the form 

Ma = -dV/da = -W'k + 2a~ 3 W 2 k - a~ 3 dW 2 /dk , 

(10) 

where we have used the form (3) for V. Notice that 
when W 2 is zero (or negligibly weak) all the hyperangular 
momenta of the form m(r a fp — rpr a ) where a ^ (3 run 
from 1 to 3N, are conserved! 

Taking the dot product of (10) with a eliminates the 
dW 2 /dk term which is purely transverse so we get the 
Virial Theorem in the form: 



2dt 2 



(Ma 2 ) = Ms 2 



dW 
da 



\w 2 =2E--^-(a 2 W) . 
a z a da 



(11) 

Multiplying by d(a 2 ) / dt / '2 and integrating Eq. (11) yields 



-Ma 2 a 2 = Ea 2 



Wa 2 - -ML 2 , 



where the final term is the constant of integration. On 
division by a 2 



M ( - C 2 \ 
Y (a 2 + -)+W(a) 



E, 



(12) 



which is the energy equation of a particle of mass M 
moving with angular momentum C and energy E in a 
hyperspherical potential W(a). Now 



— (aa + aa) = aa H —(a a) 

dt a dt 



and from (12) 



C 2 1 

a 3 M 



(13) 



(14) 



Inserting these values for a and a into equation (10) we 
obtain on simplification, multiplying by a 3 : 



Ma 2 



>da\ 



= 2W 2 k - dW 2 /dk - MC 2 k . (15) 



dt V dt J 

On writing d/dr = a 2 d/d£ this becomes an autonomous 

equation for a(r). Since (a) 2 = (aa + da) 
we may write the energy in the form 



2 - a 2 a 2 + d 2 



E 



1 



1 



= I M ( a 2 + ^+W, 
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so 

i M (^j +W2 =i MC? . (16) 

This shows that the only 'potential' in the hyper angular 
coordinates' motion is W2(a) and that the effective hyper 
angular energy in that motion is \MC? (which has the 
dimension of a 2 times an energy). In fact we showed 
in papers I and III that it was this quantity that was 
equally shared among the hyperangular momenta in the 
statistical mechanics of perpetually pulsating systems. 
(15) gives the equation of motion for a. In terms of the 
true velocities vi , V2 etc 

V^* (El, = 
dr at \ a a a J 

fa a \ 

a vi ri, v 2 r 2 ... etc , 

\ a a J 

so it is the peculiar velocities after removal of the Hubble 
flow aYi/a and after multiplication by a that constitute 
the kinetic components of the shared angular energy in 
the angular potential ^(a). 

When the m(da/dr) 2 /2 are large enough to escape the 
potential wells offered by W2/N we get an almost free 
particle angular motion of these 3N — 4 components. The 
—4 accounts for the fixing of the centre of mass and the 
removal of a from the kinetic components. If we impose 
also a prescribed total angular momentum, the number of 
independent kinetic components would reduce by a fur- 
ther three components to 3N — 7. However the peculiar 
velocities are then measured relative to a frame rotat- 
ing with angular velocity ft where ft = X~ Y • J where J 
is the angular momentum of the system. Here X is the 
inertial tensor, not to be confused with I = trace(X)/2 
introduced in Section II A 1. The constant Lagrange mul- 
tiplier is no longer f2, which is proportional to a~ 2 since 
the inertial tensor, has that dependence in pulsating equi- 
libria. The Lagrange multiplier is J = fta 2 which is 
constant during the pulsation and the peculiar velocity 

is v p i = (vi — Ui) where = ^dr^/a + J x r^/a 2 ^ see 

paper III equation (17). 1 Hereafter we specialise to the 
requirements given by Eq. (5) and (6). 

B. The 7 = 5/3 analogous fluid system 

When particles attract with both a long range force 
such as gravity or our linear law of attraction, and a 
short rage repulsion, the latter acts like the pressure of 



The above definition of u is correct. That given under equation 
(18) of paper III has — H for $1 in error as may be seen from 
equation (17). 



a fluid. Indeed short range repulsion forces only extend 
over a local region, and for large N their effect can be 
considered as pressure since only the particles close to 
any surface drawn through the configuration affect the 
exchange of momentum across that surface. In our case, 
the local r~- 3 forces come from the r~ 2 potential which 
scales as r~ 2 . A barotropic fluid with p — p 1 has an in- 
ternal energy that behaves as p 7_1 scaling like r -3 ^ 7-1 ). 
The required r~ 2 scaling gives a 7 of 5/3, i.e. a polytropic 
index of n = 3/2, the same as a non relativistic degener- 
ate white dwarf. Thus we may expect analogies between 
a 7 = 5/3 fluid with the linear long range attraction, and 
our large N particle systems. However an repulsion 
between particles is not of very short range so this fluid is 
not the same as the large TV limit of the particle system. 
The particle equilibria have the inverse cubic repulsion 
between particles balancing the long range linear attrac- 
tion, which can be exactly replaced by a linear attraction 
to the bary centre proportional to the total mass. For the 
fluid, it is the pressure gradient that balances this long 
range force. 



1. Mass profile of the static polytropic fluid 

The equilibrium of such a fluid in the presence of a 
long range force, -G*Mr, (where MG* = uo 2 ) is given 

by 

I^ = -G*Mr. (17) 
p dr 

Setting p = Kp 5 / 3 , with £ = r/r m , this yields 

P = P0 (i-ef /2 , as) 

where r is the 3D configuration space vector measured 
from the centre of mass of the system and r its modulus, 
with r m its value at the edge. 

Integrating the density gives a mass profile, M(£), 

m) = lr ( 3arcsin ^) - £( 8 £ 4 - 14 ? 2 + vVi^e) , 

(19) 

and M is the total mass. In practice for the discrete sys- 
tem of Section II A, there is a marked layering at equi- 
librium, and the "pressure" and the mass profiles depart 
from their predicted profiles as each monolayer of parti- 
cles is crossed. See Section III and Fig. (3). 

The relationship between the polytropic coefficient, ft, 
and if 2, the strength of the repulsion (entering W 2 in 
Eq. (6)) is found by identifying the internal energy of the 
corresponding fluid to the potential energy of the l/r 2 j 
coupling. In short, the former reads for a 7 = 5/3 fluid 
with density profile Eq. (17): 

V 2 = - r -^LpVWpdr = —n^rl , (20) 
Jo 7-1 128 
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while the latter reads: 

^2 



K 2 fff Tm 87rr 2 p(r)irr' 2 p(r')dpdr dr' _ 97r 4 2 4 



o 



r 2 + r' 2 — 2\irr' 



320 



(21) 

Identifying Eq. (20) and (21) gives (using M = 
7r 2 po r m/8 consistant with the density profile, Eq. (18)) 



„ = ^. 7r 4/3 M l/3 x ^ 

125 



(22) 



Similarly, requiring that V2 and V_2 balance at static 
equilibrium, where 

,2 Jl Q_2 



y_ 2 = / p(r) ^fVdr = ^PdrJL , (23) 



yields 



Po 



5G*\ 3/4 M 



, and 



, ( ^v /4 

1 5G* 



(24) 



(25) 



Eqs. (18)- (25) allow us to relate the properties of the 
crystal to the properties of the analogous fluid system. 
Notice that r m is independent of M and po is proportional 
to M. 



3. Dynamics of the rotating oscillating 7 = 5/3 fluid 

The basic pulsation of the 7 = 5/3 fluid in ro- 
tation is given by a time dependent uniform expan- 
sion/contraction plus a rotation at constant angular mo- 
mentum. Thus 



u 



-2 J x r. 



(31) 



where a(t) is the "expansion factor" of Section II A and J 
is the constant Lagrange multiplier associated with angu- 
lar momentum conservation. The acceleration involved 
in this motion are 

Du du „ du ^ ( 1 9 \ , N 

m=W +U - VU= W +V (2 U J- UX(VXU) - (32) 

Given that V x u = 2J/a 2 , putting Eq. (31) into (32) 
yields 



Du ^ far 2 1 ,~ n9 



Differentiating Eq. (8) yields: 

.. c 2 



(J 2 a . 



(33) 



(34) 



2. Figure of rotating configurations 



Equilibrium in the rotating frame with the angular 
rate, ft (not to be confused with cj, the strength of the 
harmonic potential introduced in Eq. (5)), requires that: 

-Vp = V(^ + -^ 2 ^ 2 ), where ijj = --G*Mr 2 , 
p 2 2 

(26) 

where R is the distance off axis (r 2 = R 2 + z 2 ). Since 
p = ftp 5 / 3 , it follows that, with ^ 2 = ^ 2 /(G*M), 

5«p2/3 = -^i( z 2 + (1 - tt 2 )^ 2 ) + const. (27) 



Let /) = at 2: = on axis, then 



2 E>2\3/2 

f = [i-^-(i-^A 



(28) 



Po ^0 ^0 

From this we can derive the mass, M and the moment of 
inertia, X as functions of the zo, Po 

M = Tr 2 l 2 ' and x = i MR o' ( 29 ) 

where # 2 = z 2 /(l - ^ 2 ) and z 2 = 5k/ P 2 /3 /C 2 . Eq. (28) 
generalizes Eq. (18) when the 7 = 5/3 fluid is given some 
angular momentum. The ellipticity, £, of the rotating 
configuration is given by 

e 2 = 1/(1 -S%) = 1/(1 -Q> 2 )- (30) 

Fig. (6) displays the corresponding configuration for A 7 " = 
64. 



while Euler's equation reads: 



Du 

Dt 



(35) 



Equating Eq. (35) and (33) together with Eq. (34) yields: 

V (^ ,3 + i £V -i( iXr ) 2 ) = °- (36) 

Now in the rescaled space, p = p/a 3 , r = ra and V = 
V/a, Eq. (36) reads 



0. 



(37) 



which is the equilibrium condition but with the time- 
dependently rescaled variables replacing the original 
ones. Hence it follows that in the comoving rotating 
frame, the 7 = 5/3 fluid will stratify in the same man- 
ner as Eq. (28) 2 but in the rescaled variables (Rq = 
Zq/(1 — J 2 / C 2 ) and Zq = bn/ p 2 J 3 / C 2 ). This maintains 
a constant ellipticity during the oscillation, since Ro and 
i?o are independent of a(t). 



2 which follows from integrating Eq. (37) while evaluating the in- 
tegration constant at f = 0; this integration constant could in 
principle depend on time but because J pd 3 r = M is indepen- 
dent of time, 5hip^ 3 is indeed constant 
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III. APPLICATIONS: CRYSTALLINE FORMS 
OF PULSATING EQUILIBRIA 

Let us now study the discrete form of rotating pulsat- 
ing equilibria that the systems described in Section II 
follow. Let us first ask ourselves what would the final 
state of equilibrium in which TV particles obeying Eq. (9) 
would collapse to, if one adds a small drag force in order 
to damp the motions. 



Numerical setup 

A softening scale, s, of 0.05 was used 3 so that The 
effective interaction potential reads 



K 2 



(38) 



We may escale both the repulsive, K 2 and the attractive 
strength, G* of the potential to one by choosing appro- 
priately the time units and the scale units. As a check, 
we compute the total energy of the cluster, together with 
the invariant, mC 2 /2 and check its conservation. 



A. Static Equilibrium 

The relaxation towards the equilibrium should be rela- 
tively smooth in order to allow the system to collapse into 
a state of minimal energy. We will consider here two dif- 
ferent damping forces; first (in Section III A) an isotropic 
force, proportional to —aw so that both the rotation and 
the radial oscillations are damped; or, in Section IIIB 1, 
a drag force so that the components of the velocity are 
damped until centrifugal equilibrium is reached. 



2. Larger TV limit 

Figure 3 displays both the mean mass profile and the 
corresponding section though 512 < TV < 2048 lattices; 
both the mass profile and the sections are derived while 
stacking different realisations of the lattice and binning 
the corresponding density. Note that the inner regions 
are more blurry as TV increases; indeed, the inner layers 
are frozen early on by the infall of the outer layers. 

Let us estimate the pressure profile in our crystal. In 
static equilibrium, given Eq. (17), we should have ap- 
proximately 



p(r) 



I 



uj 2 p(r)rdr 



— r^v^Ly™. (39) 



r'>r 



Eq. (39) is compared to Eq. (17) in Fig. (4). 



3. Mean density profile 

Fig. (5) shows the mean density within rescaled radius, 
£ = r/r m as a function of £ for three different simulations. 
The shell structure is very obvious but at larger TV val- 
ues it is also clear that besides the oscillations due to 
shells the mean density falls off somewhat towards the 
outside. However such a fall off is far less pronounced 
than that predicted by the polytropic model which gave 
p(£) = M(£)/(47rr^ 3 / 3 ) with M(f) given by Eq. (19). 
On further inspection, it transpired that this is due to the 
longer range of the repulsion force which is not correctly 
represented by the pressure analogy. 

Consider a continuum density distribution with a re- 
pulsive force between elements dm and dm derived from 
the potential i^2dmdm/|r — f | 2 . Then, for a spherical 
density distribution, p(f), the total potential reads 



2nr 2 p(r)dpdr 
r 2 + f 2 — 2 rrp 



2ttK 



rp(r) 



log 



dr 



1. Few particle Equilibria 

Figure 1 displays the first few static equilibria, while 
Table I lists the first 18, which includes in particular the 
regular Icosahedron for TV = 12. Strikingly, roughly be- 
yond this limit of about TV = 20, there exists more than 
one set of equilibria for a given value of TV, and the con- 
figuration of lowest energy is not necessarily the most 
symmetric. The overall structure is not far from hexag- 
onal close packed, but with a spherical layering at large 
radii. Java animations describing the crystals are found 
at http : //www. iap.fr/users/pichon/nbody.html. 



The condition of equilibrium balances against the har- 
monic attraction yields 



- — = G*Mr . 
dr 



r < r n 



(40) 



Note that all the mass, not just that inside radius r, 
contributes to the attraction for the linear law. Eq. (40) 
generates a linear integral equation for p(f). A numerical 
solution to Eq. (40) shows that the mean density distri- 
bution agrees well with the simulations, once the shell 
structure is smoothed out (see Fig. (5)). The solution to 
the integral equation is well approximated by 



P(r) = Po 1 



g^j- J , where r < r n 



(41) 



3 we also checked that our results remained the same with s = 0.01 



and beyond. The corresponding mean density within 
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Mass profile 
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FIG. 3: top: (Color online) Mean mass profile and section (top) though a N = 512, N = 1024 and N = 2048 crystal; The 
mass profile is derived while stacking about 20 realisations of the lattice and binning the corresponding density. The shell 
structure (bottom) is clearly apparent on these sections. The number of shell, N, present is consistent with the prediction of 
Section III A 2. The ordering of the inner regions is lower as N increases since the inner shells do not settle gently as they are 
disturbed by the infall of the outer layers. 



r, n, scaled to one at the centre, reads 



p(r m oe^ = 1 



3£2 ^ 

10 + 28 648 

(42) 

with £ = r/r m , and n(r) = p(r)/m, (see Fig. (5)) so that 
at £ = 1, n(l) = X = 16 651/22 680. The simulations 
have been scaled so that £ is one at the outermost parti- 
cle. That will not be at the point at which the theoretical 
smooth density falls to zero which we have called £ = 1 in 
our theoretical calculation. In practice this falls outside 
the last particle, hence we have to rescale the theory's £ 
by a factor / which is determined to make the theoretical 
mean density profile, n, fit the profile of the simulations' 
mean density. Once this scaling / has been determined 
the predicted mean density is given in terms of the ob- 
served £ by the expression Xn(/£)A/"/[4/37r(/r m ) 3 ]. 



4- Stratification in spherical shells 

In atoms, the main gradient of the potential is towards 
the nucleus, so the shell structure is dominated by the 
inner shells which have large changes in energy. In our 



systems, the global potential gradient is proportional to r 
so the largest potential differences are on the outside. As 
a consequence, the system adjusts itself so that the out- 
ermost shell is almost full, and the central region is no 
longer shell dominated. Instead of counting shells out- 
wards from the middle as in atoms, it is best to count 
shells inwards from the outside where they are best de- 
fined. For shells which are numbered from the outside, 
given that dr/n _1 / 3 (r) is the fractional increase in layer 
number, we have (using Eq. (41)) 



dr 



r . n(r) 1 / 3 



1/3 

Iq r m 
X 



1 

18 



(i - a 



(43) 

with no = po/m = 3iV/(47rr^). Conversely, the number 
of particles, N(< i) within (and including) layer i is given 

by 



1 



JLf2 + J_M _ J_,6 

10 s ' 28 s * 648 Sl 



(44) 



Solving the implicit Eq. (43) for the relative radius & of 
layer i and putting it into Eq. (44) yields the number 
of particles within layer i as a function its rank. The 



9 




Radius 



FIG. 4: (Color online) Mass profile and Pressure of an N — 
2048 simulation together with a fit of its analytic prediction 
given by Eq. (18) and (19) (see Section III A 4 and Fig. (5) 
for a discussion of the accuracy of this fit). The mass profile 
is computed via the cumulative number of particles within a 
given shell, while the pressure profile is estimated via Eq. (39). 



radii of the layers, &r m follows from inverting Eq. (43) 
for Examples of such shells are shown in Fig. (2) 
and (3), while the corresponding mass profile checked 
against Eq. (19) in Fig. (4). The number of shell present 
in Fig. (3) is consistent with the prediction of Eq. (44). 



B. Rotating Perpetually Pulsating Equilibria 

1. Rotating crystals 

The rotating crystal is achieved in steps; first, for each 
particle in the lattice, we added a velocity kick so that 
v —> v + Q x r. We then rescale the z coordinate and 
the v z (as defined along the momentum direction) of each 
particule by a constant factor. The system has then the 
ability to redistribute it momentum along the oscillating 



mean number density 




0.2 0.4 0.6 0.8 1.0 



FIG. 5: (Color online) shows the mean number density profile, 
Nn(£) for a set of about 10 N = 512, 1024 and 2048 simula- 
tions as a function of rescaled radius. The wiggles correspond 
to the shells seen in Fig. (3). The thin line corresponds to the 
prediction of Eq. (19) for the fluid system, but with an outer 
radius rescaled by the predicted radius, Eq. (25) divided by 
the measured outer radius, r m ; the dashed line corresponds 
to a uniform density solid, while the dotted dashed line cor- 
responds to the prediction of Eq. (42) (rescaled, see text). 
Clearly, this last model gives the best fit to the mean pro- 
file, which suggests that for those values of iV, the crystal 
does behave according to the equilibrium Eq. (40), and nei- 
ther like a pure solid nor as a fluid. The stars correspond to 
the predictions of the shell radii given by Eq. (43). 

particles. We then add a drag force 4 proportional to 

—a (v — X~ l • J^j , so that the components of the velocity 

are damped until centrifugal equilibrium is reached. This 
defines the rotating spheroid equilibrium. An example of 
such a configuration is shown in Fig. (6) for N = 64. The 
properties of the corresponding analogous fluid system 
are described in Section II B 2. 



2. Rotating pulsating crystals 

In order to create a pulsating configuration which pre- 
serves the shape of the rotating spheroid, we rescaled all 
the positions by some factor, A, and rescaled accordingly 
all velocities by the factor 1/A, so that angular momen- 
tum is preserved. When this special condition is met, the 



4 which can be shown to be truly dissipating energy 
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FIG. 6: (Color online) spheroidal crystal corresponding to 
the final state of a run of N — 64 with a damping force 
of the form — a (v — Half of the particles was 

painted one color, while the other half was left some other 
color. The shape preserving oscillation is achieved by rescal- 
ing all position by some factor, A, and rescaling all ve- 
locities by the factor 1/A. If the rescaling does not pre- 
serve momentum of each particle, some of the excess energy 
may go into heating the system and breaking its structure, 
which will induce mixing of the two populations (see Sec- 
tion IIIC2 and paper III). Java animations corresponding to 
the corresponding oscillating rotating crystals are found at 
http : //www . iap . f r/users/pichon/nbody . html. 
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FIG. 8: (Color online) quasi specific heat measured by in- 
crease of quasi energy with quasi temperature in units of kN; 
the quasi temperature is a 2 times the kinetic energy relative 
to the time dependent "Hubble flow" divided by 3/2 Nk; the 
quasi energy is the sum of the quasi kinetic energy plus Wi\ 
W2 is a 2 times the repulsive part of the potential energy. 
During the large amplitude pulsation of the system, the quasi 
energy is conserved. It is the quasi energy that is shared be- 
tween the different components in the statistical equilibrium. 
Since W2 is only repulsive, the system displays characteristics 
similar to a hard sphere fluid. There is a continuous transi- 
tion between the cold crystal lattice and the "free" fluid. Each 
point is derived over a set of about 20 independent runs. The 
initial condition is set up with some random motions above 
the lattice equilibrium. 



C. Thermodynamics of dissolving crystal 



1. Specific heat & evaporation 



FIG. 7: (Color online) Same as Fig. (3), but for a set of 20 ro- 
tating configuration of N = 256 particles launched according 
to the prescription given in Section IIIB1. The shell struc- 
ture is also clearly apparent on this section. The ellipticity of 
the crystal is found to be in agreement with Eq. (30). 



system oscillates and pulsates without any form of relax- 
ation (cf. Section IIIC). Note that this situation differs 
from normal modes of more classical systems, which do 
preserve the shape of a given oscillation, but might not 
involve the same particles at all times. Java animations 
describing the pulsating and rotating crystals are found 
at http : //www . iap . f r/users/pichon/nbody .html. 



How does the shell structure disappear as a function of 
temperature increase ? Fig. 8 displays the quasi-specific 
heat measured by the increase of quasi energy with quasi 
temperature in units of kN; Here the quasi temperature 
is a 2 times the kinetic energy relative to the time depen- 
dent "Hubble flow" divided by 3/2 Nk; while the quasi 
energy is the sum of the quasi kinetic energy plus W2 
given by Eq. (6). During the large amplitude pulsation 
of the system, the quasi energy is shared between the 
different components in the statistical equilibrium. Since 
W2 is only repulsive, the system displays characteristics 
similar to a hard sphere fluid. The quasi-heat capacity 
changes from the solid's 3N k to the gaseous 3/2 Nk as 
expected. 
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time time 

FIG. 9: (Color online) left panel: A?iwb/wwb, the relative 
dispersion in mixing as a function of time for a set of N = 64 
particles in a situation where the inner fraction of particles 
was rescaled along the rotation axis in order to induce some 
momentum exchange. Here three scalings (1.2, 1.6 and 2) 
were imposed (from top to bottom), corresponding to an in- 
creasing heat content which induce a faster relaxation. The 
initial condition corresponds to the prescription described in 
Section IIIB1 and Fig. (6). Right panel: Q in (t)/Q(t) and 
flout (t)/Q(t) as a function of time; again the hotter initial 
configuration (rescaled by 1.8 compared to 1.2) (full symbol) 
reaches a regime of uniform rotation more rapidly. 



2. Relaxation of the rotating pulsating configuration 

Let us split the particles within our spheroidal equi- 
librium in two sets as a function of axial radius, R, one 
corresponding to an inner ring, and one corresponding to 
an outer ring (which should initially rotate at the same 
angular rate) and rescale the z component of each inner 
particle by a factor of 2 so that it does not satisfy the 
equilibrium condition. 

A measure of the thermalisation is given by the rate 
at which the system will redistribute momentum in order 
to achieve uniform rotation again. Fig. (9) (right panel) 
represents Vt- in (t) / Ct(t) and Q out (t)/Cl(t) as a function of 
time (with Q(t) = [ft in (t) + n out (t)]/2). 



rotation axis by a factor of two. The redistribution of 
the excess energy in height will convert a fraction of it 
into heat. Fig. (9) (left panel) represents AnwB/^WB = 
RMS(n B (£) -n w (t))/RMS(n B (t) +n w (t)), with n R the 
density of WHITE particles, and ub the density of 
BLACK particles as a function of time. In practice, we 
bin the x-y plane in the comoving coordinates and esti- 
mate AnwB/^WB accordingly. 

IV. CONCLUSION 

In this paper, we have demonstrated that very special 
systems, those for which the oscillation of the system sep- 
arates off dynamically from the beheaviour of the rescaled 
variables, can form (possibly spinning) solids (albeit ones 
that pulsate in scale). We studied their phase transition 
and the corresponding "specific heat" as these structures 
melt. As expected, we found that the heat capacity halve 
to that of a free fluid, but the phase transition appears to 
occur gradually even at large N. Although we expected 
a set or regular and semi regular solids, which we did find 
for small TV, and an hexagonal close packing of most of 
the system at larger TV, we did not foresee the strong 
shell structure found. Nevertheless, we constructed a 
theory that explained this shell structure once it had 
been recognized. We also investigated the evolution of 
the lattice, as some rotation was imposed on the struc- 
ture and studied its relaxation under such circumstances, 
both in terms of mixing and for the redistribution of an- 
gular momentum. Another by-product of our study was 
that the internal energy of the corresponding barotropic 
fluid scales as p 7_1 , i.e. as (scale) -3 (^ -1 ) so that for 
7 = 5/3, it scales like (scale) -2 . Thus our investigation 
applies equally to a 7 s / 3 fluid with a quasi gravity such 
that every element of fluid attract every other in propor- 
tion to their separation. It is well known, since Newton, 
that such attractions are equivalent to every particle be- 
ing attracted in proportion to its distance to the centre of 
mass as though the total mass were concentrated there. 
With such a law of interaction, we described the possi- 
ble set of rotating pulsating equilibria it may reach, and 
found them is qualitative agreement with our simulations 
of the corresponding crystals, though the details of their 
mass profile differed: this fluid model failed to give the 
correct layering, while a more acurate model gives it cor- 
rectly. Recently [14] used simulations of such systems to 
test their smooth particle hydrodynamics codes. 



3. Mixing of the rotating pulsating configuration 

A measure of mixing is given by the rate at which the 
system becomes uniform, when it is started as two dis- 
tinct phases. Let us start by a configuration where half 
of the particles on one side of the spheroidal rotating 
equilibrium are colored in WHITE, and the other half, 
in BLACK, and rescale again the coordinate along the 
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